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Abstract —  Distributed  Radio  Frequency  (RF)  Tomography  is  a 
novel  approach  for  detecting  underground  tunnels  and  caves, 
over  relatively  wide  areas.  This  method  requires  a  set  of 
affordable  transmitters  and  receivers  randomly  deployed  above 
the  ground.  Using  the  principles  of  inverse  scattering,  it  is 
possible  to  develop  a  simplified  theory  for  imaging  below 
ground,  thus  revealing  and  locating  voids.  In  this  paper,  we 
introduce  inversion  schemes  suited  for  ground  penetrating 
radar  tomography  applied  to  sensors  that  are  randomly 
distributed  over  the  area  of  interest. 

I.  Introduction 

In  this  work,  we  address  the  incumbent  problem  of 
detection,  localization  and  monitoring  of  man-made  tunnels, 
caves,  bunkers,  weapon  caches,  underground  facilities  (UGF) 
over  relatively  wide  and  deep  regions  of  interest,  for  both 
cooperative  and  non-cooperative  (denied)  terrains.  In  military 
applications,  the  pursuit  of  “situational  awareness”  for  the 
decision  maker  is  paramount  for  successful  dominance  of  the 
scene,  especially  for  reconnaissance  /  discovery  of  UGF 
shafts  or  adits.  Similarly,  homeland  defense  and  governing 
bodies  need  different  but  equally  ubiquitous  information  to 
make  effective  decisions  regarding  disaster  response  and 
relief  activities.  Additionally,  the  engineering  and  scientific 
community  may  benefit  from  the  increased  information 
concerning  the  underground  scene,  e.g.  for  protection  of 
national  borders,  surveillance  of  sensitive  areas  (such  as 
prisons,  banks,  and  nuclear  power  plants),  mining  industry 
and  safety,  geotechnology,  environmental  engineering, 
geophysics,  archaeology  and  speleology.  According  to  the 
Layered  Sensing  philosophy  [1],  the  delivery  of  timely, 
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actionable,  trusted,  persistent  and  relevant  information  of  the 
below-ground  scene,  provided  with  affordable  architectures 
and  minimal  resource  allocation  and  human  intervention,  is 
imperative  for  the  accomplishment  of  situational  awareness. 
Although  imaging  below-ground  targets  is  presently  supplied 
by  a  plethora  of  methods,  ranging  from  seismic  to 
electromagnetic  waves,  or  from  gravity  to  optics,  from 
impedance  tomography  to  magnetotellurics,  no  technique 
complies  with  the  guidelines  prescribed  by  Layered  Sensing. 
Amid  those  techniques,  Ground  Penetrating  Radar  deserves 
particular  attention.  Recent  integration  of  diffraction 
tomography  with  classical  GPR  [2] -[8]  augmented  the 
potential  target  resolution.  However,  several  aspects  are 
challenging  the  compatibility  of  GPR  with  Layered  Sensing: 

•  Classical  GPR  resolution  is  generally  improved  by  using 
large  bandwidth,  oftentimes  requiring  high  frequencies. 
However,  due  to  soil  properties,  higher  frequencies 
experience  higher  attenuations,  and  the  increased 
frequency/bandwidth  affects  the  signal  to  noise  ratio  (SNR) 
and  intensifies  dispersion  effects. 

•  Systems  using  lower  frequencies  require  electrically  small 
wideband  radiators,  which  results  in  complex  wideband 
systems;  yet,  they  may  not  provide  adequate  resolution  to 
determine  the  geometry  of  targets. 

•  The  allocable  frequency  spectrum  is  severely  limited  by 
interference  with  unintentional  (e.g.  broadcasting  stations) 
or  intentional  (e.g.  jammers)  external  radio  sources; 
therefore,  a  reliable  system  must  operate  using  a  small, 
discrete  and  selectable  number  of  frequencies. 

•  Any  GPR  survey  requires  the  operator  to  be  above  the  area 
under  investigation. 
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•  Borehole  GPR,  which  may  increase  the  penetration  depth, 
is  generally  expensive,  subject  to  drilling  misalignments 
and,  most  important,  is  impractical  in  inaccessible 
territories. 

•  Airborne  GPR,  which  does  not  require  substantial  human 
intervention,  is  challenged  by  the  strong  reflection  due  to 
the  air-earth  interface,  which  completely  masks  the  weak 
scattered  signal  from  buried  voids. 

To  date,  no  underground  imaging  system  conforms  to 
Layered  Sensing.  Hence,  the  development  of  an  innovative 
architecture  is  mandatory. 

II.  RF  Tomography 

To  address  these  unresolved  problems,  a  new  wide  area  / 
deep  range  below-ground  imaging  of  high  contrast  dielectric  / 
conducting  anomalies,  based  upon  a  distributed  architecture, 
is  introduced  [9] -[10].  In  this  context,  the  view  (associated 
with  the  different  transmitters)  and  the  observation 
(associated  with  different  receivers)  diversity  increases  the 
information  pertaining  to  the  scene,  while  reducing  the 
necessity  of  a  large  spectral  content  of  the  probing  waveform: 
in  principle,  using  just  a  single  frequency,  range -independent 
and  sub  wavelength  resolution  images  are  feasible. 

This  distributed  scheme,  based  on  the  principles  of  RF 
Tomography ,  considers  two  separate  sets  of  N 
electromagnetic  transmitters  (Tx)  and  M  electromagnetic 
receivers  (Rx),  commonly  referred  with  the  generic  name  of 
Transponder.  These  transponders  are  placed  on  (or  in)  the 
ground  at  an  arbitrary  and  discretionary  position,  including 
random  deployment.  The  transmitter  Tx  radiates  a  known 
waveform  using  a  suitable  polarization.  The  probing  wave 
impinges  upon  a  buried  dielectric  /  conductive  anomaly,  thus 
generating  a  scattered  wave-field.  The  distributed  receivers 
collect  samples  of  the  scattered  instantaneous  electric  field, 
and  estimate  the  complex-valued  electric  field  phasor  at  their 
locations.  Subsequently,  this  information  is  relayed  to  the 
overlying  base  station.  To  ease  the  detection  process,  one 
transmitter  and  frequency  of  operation  are  activated  per  time, 
thus  simplifying  the  receiver’s  capability  to  properly  discern 
the  origin  of  the  incoming  wave-field.  For  a  given  sampling 
time,  the  used  spectrum  is  virtually  restricted  to  a  single 
frequency  component,  thus  ensuring  ultra-narrowband  system 
architecture,  low  noise  and  affordable  cost.  To  expedite  the 
set-up  and  portability  of  the  system,  sensors  are  intended  to 
be  “dart”  shaped.  Sensors  may  be  equipped  with  built-in  GPS 
for  precision  timing  and  positioning,  and  a  S-Band 
communication  link  to  transfer  the  collected  data  to  a 
overhead  base  station.  During  the  sensor  dissemination  phase, 
transponders  may  fall  in  obstructed  regions,  or  in  region 
where  clutter  is  excessive  (e.g.  Tx  and  Rx  lay  in  proximity,  or 
over  a  vegetation  layer)  and  presumably  they  may  default.  In 


any  case,  the  proposed  image  reconstruction  process  accounts 
for  the  eventual  failure  or  obstruction  of  transponders,  by 
properly  neglecting  corrupted  sensors. 

RF  Tomography  may  also  operate  using  a  discrete  set  of 
monochromatic  frequency  components,  selected  according  to 
environmental  conditions.  A  suitable  modulation  scheme  is 
the  stepped  FM;  however,  in  this  report  the  conclusions  are 
independent  from  the  type  of  modulation  implemented.  The 
choice  of  the  operating  frequencies  must  facilitate  the 
penetration  of  electromagnetic  waves  into  the  ground,  while 
simultaneously  maintaining  acceptable  detection  and 
resolution  capabilities.  Higher  frequencies  yield  better 
resolution,  but  strong  attenuation  limits  the  range  of 
operation.  Conversely,  the  corresponding  resolution 
obtainable  using  lower  frequencies  may  not  be  adequate  to 
localize  tunnels,  and  the  field  behavior  becomes  diffusive, 
thus  reducing  the  back- scattered  field.  A  suitable  range  of 
frequencies  for  this  application  is  found  to  be  the  range  1- 
15MHz.  Nevertheless,  the  final  frequency  selection  strictly 
depends  on  the  expected  target  type  and/or  depth.  Note  that  in 
HF  range  the  radiators  can  be  approximated  as  electrically 
small  dipoles  or  loops,  enabling  the  possibility  to  incorporate 
the  radiator’s  antenna  pattern  directly  in  the  signal  processing 
model.  RF  Tomography  offers  numerous  advantages  that 
technically  outperform  any  currently  available  system  for 
below-ground  imaging.  Indeed,  RF  tomography  complies  in 
full  with  the  requirements  established  by  Layered  Sensing. 
Some  reasons  are  explained  as  follows. 

Wide  area  /  Deep  Range:  the  use  of  HF  frequencies,  and 
that  the  tomographic  theory  still  holds  as  long  as  the  field 
does  not  behave  as  diffusive,  implies  that  the  probing  wave- 
front  penetrates  the  ground  at  very  large  depths.  Furthermore, 
the  burial  of  sensors,  or  their  placement  on  top  of  the  air-earth 
interface,  drastically  increases  the  amount  of  power  injected 
into  the  ground.  Assuming  ideal  conditions,  100  meters  of 
depth  and  several  thousands  of  square  meters  of  area  can  be 
theoretically  surveyed.  Due  to  the  natural  spread  of  sensors 
over  a  wide  area,  the  system  is  well  suited  for  horizontal 
prospecting.  The  classical  GPR  diffraction  tomography  was 
unable  to  include  depolarization  effects  that  normally  occur 
when  the  anomaly  is  considered  high-contrast.  Furthermore, 
diffraction  tomography  is  not  applicable  when  the  host 
medium  changes  the  polarization  of  the  probing  wave,  as  in 
the  case  of  the  air-earth  interface.  With  the  introduction  of  the 
generalized  diffraction  tomography ,  proposed  in  [7]-[15],  the 
depolarization  effect  is  accounted,  thus  yielding  better 
reconstructed  images.  Moreover,  the  proposed  forward  model 
formulation  considers  the  possibility  of  the  sensors  to  be  in 
the  near  field,  thus  allowing  larger  investigation  areas,  from 
very  shallow  to  very  deep  regions.  Therefore,  the  typical 
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problem  of  “blind”  region  that  affects  any  radar -based  system 
is  virtually  eliminated. 

Actionable  Deliverables:  the  forward  model  for  RF 
tomography  and  the  proposed  generalized  diffraction 
tomography  have  been  developed  by  neglecting  part  of  the 
information  concerning  the  value  of  the  dielectric  permittivity 
and  the  shape  of  the  targets,  while  preserving  the  information 
concerning  the  location  and  the  overall  extension  of  the 
object.  Therefore,  the  proposed  tomographic  reconstruction 
procedure  has  been  tailored  to  ease  the  process  of  detection , 
identification,  location,  trace  and  reconnaissance  of 
underground  bodies,  thus  providing  the  decision  maker  only 
with  intelligible ,  relevant  and  actionable  information. 

Minimal  Human  Intervention:  RF  tomography  is 
particularly  suited  for  remote  probing  without  the 
involvement  of  human  resources.  Considering  a  dangerous 
area,  where  humans  would  be  in  danger,  one  can  still 
surround  the  area  of  interest  with  transponders  (for  example 
along  a  circle)  and  RF  tomography  can  still  provide  relevant 
information  of  below  ground  entities.  Furthermore,  in  the 
event  that  the  area  is  completely  inaccessible,  sensors  can  be 
disseminated  from  an  airborne  vehicle,  and  below-ground 
images  can  be  reconstructed  in  a  very  short  amount  of  time. 

Real  Time:  the  time  required  by  RF  tomography  to 
reconstruct  an  actionable  image  of  below  ground  targets  is  the 
mere  summation  of  the  time  needed  for  deploying  sensors 
above  the  region  of  interest,  and  the  time  requested  for  data 
collection  and  retrieval.  If  sensors  are  dropped  from  an 
aircraft,  and  the  data  processing  is  performed  via  a  high 
performance  computer,  then  reconstructed  images  can  be 
available  in  terms  of  minutes,  and  revisit  time  is  measurable 
in  terms  of  seconds. 

Affordability :  the  design,  manufacturing,  testing  and 
realization  of  transponders  in  the  HF  band  are  relatively 
simple  and  inexpensive  procedures,  and  precise  location, 
timing,  and  phase  coherence  can  be  guaranteed  by  the  well 
established,  standardized,  and  low-cost  GPS  technology. 

Agility:  differently  from  common  geophysical  survey 
methods,  where  regular  grids  are  necessary  to  sample  the 
data,  RF  tomography  is  developed  in  a  way  that  the  location 
of  sensors  is  discretionary,  and  it  may  be  completely  random. 
Randomness  of  the  sensor  deployment  does  not  preclude  the 
system  operability:  at  most,  if  sensors  are  poorly  arranged  in 
the  ground,  a  slight  degradation  of  the  reconstructed  image 
may  be  expected.  By  accounting  for  the  possibility  of  default 
of  transponders,  the  tomographic  inversion  procedure  can  be 
considered  modular,  and  the  sensing  can  be  persistent. 
Furthermore,  RF  tomography  is  suited  for  detection  of  targets 
of  any  shape  and  from  any  sensor  distribution.  In  fact, 
conventional  radar  systems  are  not  able  to  detect  targets 


whose  backscattered  RCS  is  equal  to  zero.  In  order  to 
increase  the  detection  probability,  different  waveforms, 
multi  static  solutions  or  auto -focusing  need  to  be  implemented 
in  the  system.  In  RF  Tomography,  any  particular  value  of  the 
RCS  is  considered  useful  information,  including  the  case 
when  RCS  =  0 .  Therefore,  RF  Tomography  is  virtually  able 
to  detect  targets  that,  from  some  particular  points  of  view, 
may  have  zero  backscattered  field,  and  use  this  information  to 
accurately  determine  properties  of  the  target. 

Spectral  Dominance:  theoretically,  RF  tomography 
features  the  highest  number  of  diversities  among  any  other 
void  detection  system.  In  fact,  RF  tomography  comprises 
frequency,  polarization,  waveform,  code,  modulation, 
antenna  pattern,  view,  and  observation  diversities.  RF 
tomography  works  in  principle  with  a  monochromatic 
waveform.  Indeed,  with  the  use  of  generalized  diffraction 
tomography,  the  more  the  system  is  narrowband,  the  better 
the  performance  is.  Since  RF  Tomography  operates  at  a 
single  frequency,  it  is  not  affected  by  frequency  dispersion. 

Resolution:  The  use  of  generalized  diffraction  tomography 
in  our  proposed  system  guarantees  the  sub -wavelength  and 
range  independent  resolution.  Current  research  concerning 
compressive  sensing,  time  reversal  and  MUSIC  algorithms 
may  further  improve  the  capability  of  discrimination  among 
closely  spaced  targets,  thus  achieving  the  so-called  super¬ 
resolution. 

Robustness  to  Noise:  in  RF  tomography  the  concept  of 
“sample”  is  more  elaborated  than  in  classical  radar  theory, 
where  samples  are  collected  at  different  times,  and  then 
averaged,  to  increase  the  SNR.  In  our  case,  a  “sample”  may 
be  defined  as  the  particular  combination  of  transmitter, 
receiver,  frequency  (or  waveform),  polarization,  and  time. 
Following  this  definition,  for  a  given  time,  frequency  of 
operation  and  polarization,  we  may  expect  a  SNR  increase 
proportional  to  the  square  root  of  the  total  number  of  Tx  and 
Rx.  By  time  averaging  these  results,  the  expected  SNR  can  be 
arbitrarily  increased  according  to  the  operator’s  demand. 
Although  this  unbounded  improvement  is  not  reached  in  real 
cases  (because  clutter  cannot  be  modeled  as  Gaussian  noise), 
it  still  gives  an  insight  of  the  tremendous  robustness  to 
noise/clutter  of  RF  tomography. 

Persistent  Sensing:  the  peculiar  distributed  sensing  in  RF 
tomography  provides  the  user  with  profound  situational 
awareness  of  the  underground  scene.  Particularly,  the 
persistence  is  guaranteed  by  the  multitude  of  sensors,  and  the 
failure  of  some  Tx/Rx  units  does  not  preclude  the  overall 
image  reconstruction  process. 
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Figure  1 .  Geometric  representation  of  RF  Tomography 


III.  Forward  Model 

To  derive  a  general  framework  for  RF  Tomography 
inversion  algorithms,  a  suitable  model  for  the  electric  field 
(closely  representing  the  actual  scattering  process)  needs  to 
be  defined.  The  3D  geometry  of  the  problem  is  depicted  in 
Fig.  1.  For  simplicity,  a  single  operating  frequency  /  is 

adopted,  but  extension  to  the  multi -frequency  operation  is 
straightforward.  Under  the  monochromatic  assumption,  the 
ground  is  modeled  as  a  semi-space  homogeneous  medium 

with  relative  dielectric  permittivity  SD  ,  conductivity  <J D  , 

and  magnetic  permeability  //0  .  The  targets  (i.e.  tunnels  or 

voids)  are  assumed  to  reside  inside  the  investigation  domain 
D,  which  is  fully  enclosed  in  the  ground.  The  sources  are  N 

electrically  small  dipoles  (of  length  A  f )  or  loops  (of  area 

Ar )  fed  with  current  V ,  and  located  at  position  Yln  (view 

diversity),  and  having  (electric  or  magnetic)  dipole  moment 
direction  and  polarization  expressed  by  the  complex  valued 

unit- norm  vector  .  For  each  transmitting  antenna,  the 

scattered  field  is  collected  by  M  receivers  (observation 
diversity),  located  at  points  in  space,  and  having  moment 

direction  and  polarization  .  The  unknowns  of  the 
problem  are  the  relative  dielectric  permittivity  profile 
£r(r')  and  the  conductivity  profile  cr(r')  inside  the 

investigation  domain  D.  The  resulting  inverse  problem  is  cast 
in  terms  of  the  unknown  permittivity  contrast  function: 


The  total  instantaneous  electric  field  experienced  at  the 
receiver  in  time  domain  has  been  derived  in  [10],  [16] -[17]: 

'■'(CCO'Keje  '“[C*o  * 

J1J0:  '  §  (r; , r  •)]  •  [g  (r  \  r; )  •  a;  ]  (r  •)  dr '  (2) 

D 

+fia:-G(r:,r;)X  +  //(r:,r; )]}  +  «(,) 

The  term  strm  •  G  (r^ , •  2tn  in  eq.  (2)  represents  the  direct 

coupling  between  Tx  and  Rx,  and  it  can  be  considered  a 

source  of  deterministic  clutter.  The  term  H  (r^r^ 

represents  the  non-linear  contribution  due  to  higher-order 
Born  series  terms.  In  the  reconstruction  process, 

H  (r^r^  is  considered  as  unknown  static  bias  (clutter)  to 

the  desired  signal.  The  random  variable  N(t}  represents  the 

thermal  noise  and  the  quantization  errors,  and  can  be  treated 
as  a  Gaussian  process  with  zero  mean. 

RF  Tomography  presents  the  challenging  problem  of 
mitigation  of  the  direct-path  coupling  between  transponders. 
A  simple  strategy  to  mitigate  the  direct-path  coupling  consists 
in  the  constrained  minimization  of  the  electric  field  radiated 
from  the  transmitting  antenna  measured  at  the  position  of  a 
particular  receiver.  Mathematically,  the  constrained 
minimization  problem  is  recast  in  the  following  form:  for 
each  n,m  pair, 


minimize:  ^(Cr^-a^ 
subject  to:  (a' )T  -a'  =1 


(3) 


This  constrained  minimization  problem  can  be  easily  solved 
by  constructing  a  series  of  Lagrangian  forms: 


A. 


-A 


W' 


•a,  -1 


(4) 


where  G  *  denotes  the  Hermitian  of  G .  By  imposing 
VA  =  0  we  obtain: 


es{r')  =  er 


q~(r’)-g'o 

27Tfs0 


q*(<yn)-q(<yn)<=K  (5) 

(1) 

Therefore,  the  an  direction  that  minimized  the  power  at  a 
desired  location  is  the  eigenvector  associated  with  the 
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smallest  eigenvalue  of  the  (positive  semidefinite)  matrix 

G*G  .  Similarly,  we  can  apply  this  minimization  technique 
at  the  receiver  side:  by  defining  the  vector  quantity 

Emm  =§(r>')-a'  (6) 


as  the  electric  field  obtained  when  a  tn  is  chosen  according  to 
(3),  a  second  minimization  problem  can  be  formulated: 


minimize 


K) 

subject  to:  (a,rn)T  < 


(7) 


The  minimization  is  achieved  when  a rm  is  chosen  to  be  the 
eigenvector  corresponding  to  the  smallest  eigenvalue  of  the 
matrix  Emin  ■  E*min . 


As  tested  via  numerical  analysis,  the  application  of  these 
strategies  guarantees  an  acceptable  minimization  of  the 
received  clutter,  thus  the  total  electric  field  (in  phasor  form) 
can  be  approximated  with  the  scattered  contribution  from 
targets  inside  the  region  D, : 


E  (< ,  r; )  =  Es  (r' ,  r’m )  =  L  (et  (r '))  =  Qkl  x 
I':  g(r:,r')]{§(r',r;)a;]£,(r')*’  <8) 

D 

IV.  Inversion  Strategies 

The  inversion  algorithms  are  valid  assuming  that  the  direct 
path  contribution  and  Gaussian  noise  are  substantially 
mitigated,  so  that  their  contribution  can  be  considered  as 
perturbations  on  the  measured  data;  therefore,  the  field 
received  by  the  sensors  is  assumed  to  be  the  one  described  by 
the  forward  model  in  eq.  (8).  A  way  to  compute  the  contrast 
function  using  (8)  is  to  perform  a  numerical  inversion  of  the 
integral  operator.  The  sampled  field  can  be  collected  in  a 

vector  E  =  ji? [y^  ,  r^J  ,  of  length  NM  ,  and  the  region  D 

can  be  discretized  in  K  voxels,  each  one  located  at  position 
Yk  ’ :  the  contrast  function  can  be  embodied  in  a  column 

vector  &  =  {s5  (rfc  of  length  K.  After  this  discretization, 
eq.  (8)  can  be  rewritten  in  matrix  form  as: 


E  = 


(9) 


conditioned.  The  operator  L  is  commonly  reaching  values  of 

condition  number  K  greater  than  106  .  This  leads  to  artifacts 
in  the  reconstruction  process,  particularly  exacerbated  when 
noise  (thermal,  external,  quantization),  clutter  or  coupling  is 
impinging  upon  the  receivers. 

According  to  the  accuracy  required  from  the  system,  several 
inversion  strategies  are  proposed  [10],  [16] -[17]: 

•  Levenberg-Marquardt  (LM)  regularization  procedure. 
This  method  is  relatively  accurate  for  any 
environmental  condition  and  it  is  robust  in  presence 
of  noise.  It  requires  a  proper  choice  of  a 
regularization  parameter. 

•  Truncated  Singular  Value  Decomposition  (TSVD). 
This  method  is  also  relatively  accurate  in  any  scenario 
and  moderately  resistant  to  noise  interference.  This 
method  also  offers  deeper  insight  into  the  physics 
behind  the  reconstruction,  and  the  output  can  be 
easily  adjusted  by  properly  selecting  the  number  of 
meaningful  singular  values.  The  number  of  retained 
singular  values  in  TSVD  plays  the  same  role  of  the 
LM  regularization  parameter. 

•  Back  Propagation  approach.  This  method  works 
properly  only  when  the  operator  L  is  well 
conditioned.  This  implies  that  it  can  be  used  only  for 
particular  configurations  and  when  the  SNR  is 
relatively  high.  However,  the  computational  time  is 
drastically  reduced. 

•  Weighted  LM  and  Weighted  TSVD.  These  strategies 
are  natural  extensions  of  LM  and  TSVD  methods  that 
include  a  proper  weighting  function  of  the  sensors,  in 
order  to  equalize  the  contribution  due  to  a  particular 
Tx  and  Rx  pair.  If  speed  is  not  a  concern,  weighted 
LM  or  weighted  TSVD  are  the  superlative  choices. 

•  Fourier-Bojarski  approach.  This  is  the  fastest 
inversion  scheme  and  it  is  suited  for  far-field  probing 
and  quasi-lossless  soils.  It  privileges  speed  instead  of 
accuracy. 

•  Sub-Space  Models ,  suited  for  increasing  the 
resolution  when  very  low  frequencies  are  used  and 
the  targets  can  be  approximated  as  being  point-like. 

•  Compressive  Sensing,  suited  when  the  number  of 
transmitters  and  receivers  (i.e.  the  number  of 
samples)  is  limited. 

The  mathematical  description  of  the  first  four  methods  is 
presented  in  [10],  [16-17].  In  this  work,  we  will  discuss  more 
in  detail  the  last  3  approaches. 


Hence,  the  problem  is  to  invert  the  relation  (9).  In  eq.  (9),  L 
is  a  matrix  with  dimensions  NM  x  K ,  which  is  typically  ill 
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A.  Weighted  Regularization  Methods 

Weighted  methods  are  based  upon  the  assumption  that  the 
insertion  of  weighing  matrices  opportunely  defined  may 
equalize  the  information  content  provided  by  each  Tx-Rx 
pair.  For  instance,  the  weighted  version  of  the  LM  method 
can  be  reformulated  as  follows  [16-17]: 


The  first  step  is  to  use  all  transmitters  simultaneously.  We  fill 
the  multistatic  response  matrix  in  the  following  manner: 


K(r\)  =  klQ 


Gtx\Grx\ 


GtxnGrx] 


Gtx\GrxM 


GrxnGrxM 


,  (15) 


&(/?)  =  (l’W^L  +  /?W; )-‘ (L*W \E  +  J3 )  (10)  where 


where  &  represents  the  known  dielectric  anomalies 
embedded  in  region  D  ,  and 


Gr„  =G(r 
C,_  =G(r',r',) 


(16) 


W£=diag(LX)1/2 
W£  =  diag  (LL*  )‘/2 


(ID 


We  define  “Green’s  function  vector”  (or  steering  vector)  the 
vector  entries  in  (15): 


(12) 


s{Ck)  =  [GTx,GKx]  G,xNGrJ  (17) 


are  (diagonal)  the  weighting  matrices,  opportunely  defined  in 
order  to  minimize  the  sensitivity  of  the  system  [2]. 

In  most  cases,  the  weighting  matrices  have  a  small  dynamic 

range.  Therefore,  we  can  approximate  W£  =  al  in  the  matrix 
inversion  term  of  eq.  (10).  If  we  perform  the  singular  value 
decomposition  [11]  of  a  weighed  matrix  Lw  defined  as: 


= wel = usv* 


(13) 


we  can  rewrite  (10)  in  terms  of  (13),  obtaining: 

4  =  (l;l, + /?w;  )_1  (l;  w ee + /?w;£;’ ) 

=  (vs*sv*  +a2p\y  (YS*U*W££  +  y9W^0)  (14) 


:  Vdiag 


1 


s?+a2/3 


(SX'W^  +  V^W^0) 


where  s.  represent  the  i-th  singular  value  of  Lw .  The 
advantage  obtained  in  applying  eq.  (14)  is  that  the  contrast 
function  as  function  of  [3  can  be  computed  via  (fast)  matrix 
multiplications. 

B.  Time  Reversal  and  Sub-Space  Method 
Sub- space  methods,  such  as  the  MUSIC  algorithm,  are 
appealing  when  the  probing  frequency  is  very  low,  where 
improving  resolution  is  the  main  challenge.  The  lower  the 
frequency,  the  more  the  tunnels  resemble  point  scatterers, 
therefore  imaging  via  the  MUSIC  algorithm  may  be  possible. 


It  has  been  shown  [18]  that  the  eigenvectors  of  the  matrix 

T  =  KhK  (18) 


corresponds  in  a  one-to-one  manner  to  different  targets  (when 
targets  are  well-resolved).  In  particular,  the  wavefront 
generated  by  the  array,  when  excited  by  one  of  these 
eigenvectors,  focuses  on  the  associated  target  so  that  a 
synthetic  image  of  the  target  locations  is  easily  computed. 
This  procedure  is  commonly  referred  to  as  “ Time  Reversal ” 
[18].  However,  it  is  known  that  time  reversal  shows  poor 
performance  when  the  array  is  sparse  (which  is, 
unfortunately,  our  case).  A  way  to  improve  the  reconstruction 
quality  and  lead  to  super  resolved  images  is  to  apply  the 
MUSIC  algorithm  to  the  matrix  T.  The  MUSIC  algorithm 
makes  use  of  the  fact  that  the  time -reversal  matrix  T  is  a 
projection  operator  onto  the  subspace  spanned  by  the 
complex  conjugates  of  the  Green  function  vectors  of  K  (the 
signal  sub-space)  and  that  the  noise  subspace  is  spanned  by 
the  eigenvectors  of  T  having  zero  eigenvalue.  It  then  follows 
that  the  complex  conjugate  of  each  Green  function  vector 
must  be  orthogonal  to  the  noise  subspace  and,  in  particular,  to 
the  eigenvectors  of  the  T  matrix  [18].  We  define  the  MUSIC 
pseudo-spectrum  as: 


1 

Z|(w(r'd)|2 


(19) 


where  are  the  eigenvectors  of  T  having  zero  eigenvalue. 

According  to  MUSIC  algorithm,  peaks  of  D  reveal  the  actual 
target  locations. 
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C.  Compressive  Sensing 

Compressive  sensing  (CS)  is  a  concept  formally  introduced 
by  Candes  et  al.  [19],  and  Donoho  [20]  to  overcome  the 
Shannon/Nyquist  restriction  on  the  sampling  period. 
Recently,  the  scientific  community  has  discovered  its 
potentialities  to  many  applications,  including  imaging.  In 
fact,  problems  in  compressive  sensing  are  generally 
formulated  as  constrained  minimization  of  matrix  equations. 
Using  our  notation,  a  typical  compressive  sensing  problem  is 
recast  in  terms  of  the  following  minimization: 


minimize  \\sA 

ll^  lip 

subject  to  Lsj  =  Es 


(20) 
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Figure  2.  Sensor  disposition  for  the  case  1 . 


where  the  norm  used  for  the  minimization  is  represented  by  p. 
With  respect  to  others  “standard”  implementation  of  CS,  we 
are  not  building  a  sampling  matrix,  but  we  are  using  the 
matrix  L  as  it  is  obtained  from  the  forward  model.  When  p=2 
the  minimization  problem  is  identical  to  the  classical 
Tikhonov  regularization  procedure.  However,  when  p=l  or 
p—0  the  solution  of  the  constrained  minimization  shows  the 
peculiar  feature  of  having  a  greater  number  of  zero -valued 
pixels.  This  leads  to  high-contrast,  sharper  and  low-blurred 
images  with  respect  to  classical  regularization  procedures. 
This  result  can  be  used  also  in  a  different  way:  using 
compressive  sensing,  the  same  image  quality  is  expected 
using  fewer  samples,  which  in  turns  may  save  the  system 
complexity  and  cost.  Three  algorithms  typical  of  compressive 
sensing  are  investigated:  Basis  Pursuit  [21]-[22],  Orthogonal 
Matching  Pursuit  [23],  and  Homotopy/LASSO  [24]-[25]. 

V.  Simulations 

Simulations  are  performed  using  two  different  environments: 

7 )  Random  distribution.  The  simulation  intends  to  mimic 
the  case  in  which  sensors  are  randomly  deployed  above  an 
area  of  interest.  Targets  are  expected  to  be  located  below  the 
transponders. 

2)  Close-In  Sensing.  The  simulation  intends  to  replicate 
the  case  in  which  a  denied  area  needs  to  be  explored. 
Therefore,  sensors  are  surrounding  the  area  of  interest,  and 
images  are  reconstructed  by  exploitation  of  lateral  waves. 
Preliminary  results  are  shown  in  Figs.  4-5,  where  the  two 
tunnels  can  be  easily  recognized.  Details  regarding  the 
simulation  are  found  in  [10].  Further  results  and  comparisons 
will  be  shown  at  the  time  of  the  conference. 


X  Range  [m] 


Figure  4:  Reconstructed  image  using  TSVD  and  scattered  field  data 
generated  using  the  forward  model  in  Section  II  and  the  geometry  depicted  in 

Fig -2 
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